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Abstract. This paper presents some results on a well-known problem in Algebraic Sig- 
nal Sampling and in other areas of applied mathematics: reconstruction of piecewise-smooth 
functions from their integral measurements (like moments, Fourier coefficients, Radon trans- 
form, etc.). Our results concern reconstruction (from the moments or Fourier coefficients) 
of signals in two specific classes; linear combinations of shifts of a given function, and 
"piecewise ZJ-finite functions" which satisfy on each continuity interval a linear differential 
equation with polynomial coefficients. In each case the problem is reduced to a solution 
of a certain type of non-linear algebraic system of equations ("Prony-type system"). We 
recall some known methods for explicitly solving such systems in one variable, and provide 
extensions to some multi-dimensional cases. Finally, we investigate the local stability of 
solving the Prony-type systems. 

1. Introduction 

It is well known that the error in the best approximation of a C'^-function / by an A^- 
th degree Fourier polynomial is of order The same holds for algebraic polynomial 

approximation and for other basic approximation tools (see e.g. [21, Chapters IV, VI], 
and [221 Vol.1, Chapter 3, Theorem 13.6]). However, for / with singularities, in particular, 
with discontinuities, the error is much larger: its order is Considering the so-called 

Kolmogorow A^-width of families of signals with moving discontinuities one can show that 
any linear approximation method provides the same order of error, if we do not fix a priori the 
discontinuities' position (see [H], Theorem 2.10). Another manifestation of the same problem 
is the "Gihhs effect" - a relatively strong oscillation of the approximating function near the 
discontinuities. Practically important signals usually do have discontinuities, so the above 
feature of linear representation methods presents a serious problem in signal reconstruction. 
In particular, it visibly appears near the edges of images compressed by JPEG, as well as in 
the noise and low resolution of the CT and MRI images. 

Recent non-linear reconstruction methods, in particular. Compressed Sensing ([HI [7]) and 
Algebraic Sampling (P|19l[23l[T0l[T5]), address this problem in many cases. Both approaches 
appeal to an a priori information on the character of the signals to be reconstructed, assum- 
ing their "simplicity" in one or another sense. Compressed sensing assumes only a sparse 
representation in a certain (wavelets) basis, and thus it presents a rather general and "uni- 
versal" approach. Algebraic Sampling usually requires more specific a priori assumptions 
on the structure of the signals, but it promises a better reconstruction accuracy. In fact, 
we believe that ultimately the Algebraic Sampling approach has a potential to reconstruct 
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"simple signals with singularities" as good as smooth ones. The most difficult problem in this 
approach seems to be estimating the accuracy of solution of the nonlinear systems arising. 

Our purpose in this paper is to further substantiate the Algebraic Sampling approach. On 
one hand, we present two algebraic reconstruction methods for generic classes of signals. The 
ffist one is reconstruction of combinations of shifts of a given function from the moments 
and Fourier coefficients (Section [2]). The second one concerns piecewise D- finite moment 
inversion (Section H]). On the other hand, we consider some typical nonlinear systems arising 
in these reconstruction schemes. We describe the methods of their solution (Section |3]), and 
provide some explicit bounds on their local stability (Section [5]). We also present results of 
some numerical experiments in Section |6l 

Our ultimate goal may be stated in terms of the following conjecture (which seems to be 
supported also by the results of |9l UHl HH [301 123]): 

There is a non-linear algebraic procedure reconstructing any signal in a class of piecewise 

-functions (of one or several variables) from its first N Fourier coefficients, with the 
overall accuracy of order This includes the discontinuities' positions, as well as the 
smooth pieces over the continuity domains. 

Recently [1] we have shown that at least half the conjectured accuracy can be achieved. 
However, the question of maximal possible accuracy remains open. Our results presented in 
this paper can be considered as an additional step in this direction. 

2. Linear combinations of shifts of a given function 

Reconstruction of this class of signals from sampling has been described in [HI [19]. We 
study a rather similar problem of reconstruction from the moments. Our method is based on 
the following approach: we construct convolution kernels dual to the monomials. Applying 
these kernels, we get a Prony-type system of equations on the shifts and amplitudes. 

Let us restate a general reconstruction problem, as it appears in our specific setting. We 
want to reconstruct signals of the form 

N 

(1) F(x) = ^a,/(x-x^) 

i=i 

where / is a known function of x = (xi, . . . , Xa) G M"^, while the number of the shifts and 
the form ([1]) of the signal are known a priori. The parameters aj, x^ = {x{, . . . j = 

1,. . .N are to be found from a finite number of "measurements", i.e. of linear (usually 
integral) functionals like polynomial moments, Fourier moments, shifted kernels, evaluation 
over some grid etc. 

In this paper we consider only linear combinations of shifts of one known function /. 
Reconstruction of shifts of several functions based on "decoupling" via sampling at zeroes 
of the Fourier transforms of the shifted functions is presented in |22]- 

In what follows x = (xi, . . . ,Xd),t = (ti, . . . ,td) G M'^, j is a scalar index, while k = 
{ki, . . . ,kd), i = {ii,...,id) and n = {rii, . . . ,nd) are multi-indices. Partial ordering of 
multi- indices is given hj k < k' kp < k'^, p = 1, . . . , d. 

Assume that a multi- sequence of functions ip = {(/^^(t)}, t G M'^, > (0, . . . , 0) is fixed. 
We consider the measurements Hk{F) provided by 
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(2) /ifc(F) = J F{t)Mt) dt, A; > (0, . . . , 0). 

Our approach now works as follows: given / and f = {^k{t)} we now try to find an "/- 
convolution dual" system of functions i/j = ipk{t) in a form of certain "triangular" linear 
combinations 



(3) ipk{t) = J2 C*,fcV^*W, A;> (0,...,0). 

0<i<fc 

More accurately, we try to find the coefficients Ci^k in © in such a way that 

(4) [ f{t-x)Mt)dt = ifk{x). 



We shall call a sequence ifj = {'ipkit)} satisfying ([3j), (j4]) / - convolution dual to ip. Below 
we shall find explicitly convolution dual systems to the usual and exponential monomials. 

We consider a general problem of finding convolution dual sequences to a given sequence of 
measurements as an important step in the reconstruction problem. Notice that it can be gen- 
eralized by dropping the requirement of a specific representation ([3]): ipk{t) = Yl'i=o Ci^k'^i{t). 
Instead we can require only that J f{t)ijjk{t) dt be expressible in terms of the measurements 
sequence jik- Also ipk in (jl]) can be replaced by another a priori chosen sequence ijk- This 
problem leads, in particular, to certain functional equations, satisfied by polynomials and 
exponential functions, as well as exponential polynomials and some kinds of elliptic functions 
(see [28j). 

Now we have the following result: 

Theorem 1. Let a sequence -0 = ^|Jk{t) be f -convolution dual to ip. Define Mk by Mk = 
J2o<i<k Ci,kfJ'i- Then the parameters aj and in (QP satisfy the following system of equations 
( "generalized Prony system"): 

N 

(5) Yajipk{x^) = Mk, A; > (0,...,0). 

Proof. By definition of Mk and via ([3]) and (j4]) we have for each A; > (0, . . . , 0) 



Mk=Y^ Ci,kfii= I F{t) J2 Ci.mWdt 

0<j<A; 0<i<fc 
» JV „ N 

/ Fit)Mt)dt = J2<^j / fit-x^)Mt)dt = Y, 

.7 = 1 .7 = 1 



aj(pk{x^] 



This completes the proof. □ 

According to Theorem [H in order to reduce the reconstruction problem with the mea- 
surements ([2]) and for signals of the a priori known form ([T]) to a solution of the generalized 
Prony system we have to find an / - convolution dual system ip = ipkit) to the measurements 
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kernels ^. In fact we need only the coefficients Cj^fc. Having these coefficients, we compute 
Mfc and get system ([5]). 

Solvability of system ([5]) and robustness of its solution depend on the measurements kernels 
ip. Specific examples are considered below. 

2.1. Reconstruction from moments. We are given a finite number of moments of a signal 
-^(^) = ^Y^]=\ (^jf{^ ~ ^'^) in ([1]) in the form 



(6) nik 



j F{t)t^dt, k=ih,...,ka) > (0,...,0). 



So here (^fe(x) = x\ 



■ x/ for each multi-index k. We look for the dual functions ^|Jk 



satisfying the convolution equation 



(7) 



fit- x)Mt)(it = x'' 



for each multi-index k. To solve this equation we apply Fourier transform to both sides of 
d?]). Assuming that /(O) 7^ and that /(w) has the derivatives of all the orders at we find 
(see [25) that there is a unique solution to ([7]) provided by 



(8) 

where 



k,l 



2ny 



Kk 



Qk—l 


1 




a;=0 /(W)_ 



The assumption /(O) 7^ is essential in the construction of the /-convolution dual system 
for the monomials as well as for other measurement kernels (as well as in the study of the 
shifts of a given function in general). The above calculation can be applied, with proper 
modifications, in more general situations (see [28]). On the other hand, the assumption of 
differentiability of / at zero is not very restrictive, in particular, if we work with signals with 
finite support. 



Returning to the moments reconstruction, we set the generalized polynomial moments Mk 



as 
(9) 



M, 



kim 



Kk 



and obtain, as in Theorem [H the following system of equations: 



(10) 



N N 



^aj(x-') 



Mk, A; > (0,...,0). 



This system is called "multidimensional Prony system". It appears in numerous problems 
of theoretical and applied mathematics. In Section |3] below we recall one of the classical 
methods for its solution in one- dimensional case, and describe, under certain assumptions, 
a method for its solution in several dimensions. Local stability of the solutions of the one- 
dimensional Prony system is discussed in Section O 

4 



2.2. Fourier case. The signal lias the same form as in section I^TTl F{x) = ~ 
x^). The measurement kernels we now choose as the harmonics ^k{x) = e*^^'. So our 
measurements are the Fourier coefficients Ck{F) = f{k) = J F{x)e^''^ dx. 

Let us assume now that / satisfies the condition f{k) ^ for all integer k. Then the Fourier 
harmonics turn out to form essentially /-self-dual system: indeed, taking '^k{x) = j^e*'^^ 

we get immediately that J f{t — x)il>kit) dt = J f{t — x)^j--e^^^ dt = 77^6*^^' = (fkix). 
According to our general scheme we put now Mk = -i77TCk{F) and get a system of the form 

^ N N 

(11) c,(F) = M, = 5^a,e^'^-^ =^a,.(^y, fc>(0,...,0), 
J\i^) j=i j=i 

where p' = e^^\ This is once more a multidimensional Prony system as in ( ITOl) . with the 
nodes on the complex unit circle. 

2.3. Further extensions. The approach above can be extended in the following directions: 

(1) Reconstruction of signals built from several functions or with the addition of dilations 
also can be investigated (a perturbation approach where the dilations are approxi- 
mately 1 is studied in [26]). 

(2) Further study of "convolution duality" in [26] provides a certain extension of the class 
of signals and measurements allowing for a closed - form reconstruction. 

3. Prony system 

3.1. One-dimensional case. Let us start with the classical case of one-dimensional Prony 
system: 

N 

(12) ^aj{x^f = ruk, x^ eR, A; = 0, 1, ... . 

i=i 

This system appears in many branches of mathematics (see [201) • There are several solution 
methods available, for instance direct nonlinear minimization (see e.g. [H]), the polynomial 
reahzation (original Prony method, [21]) or the state-space approach ([25]). Let us give a 
sketch of a method based on Pade approximation techniques which is rather close to the 
original Prony method. To simplify the presentation we shall assume that we know a priori 
that all the nodes x^ are pairwise distinct. General case is treated similarly ([22], see also 

Consider a "moments generating function" I{z) = YlT=o''^kZ^ , z G C. Assuming the 
equations ffT^ are satisfied for all = 0, 1, . . . and summing up the geometric progressions 

we get 



(13) 



N 



Hz) 
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X^ z 



So I{z) is a rational function of degree tending to zero as z ^ oo. The unknowns aj and 
i in f|T2l) are the poles and the residues of I{z), respectively. 

Now in order to find I{z) explicitly from the first 2N moments mo, mi, . . . , "we use 
the Pade approximation approach (see [22]): write I{z) as with polynomials P{z) = 

Ao + Aiz H h Am-iZ^~^ and Q{z) = Bq + Biz H h -Batz^ of de grees — 1 and N , 

respectively. 

Multiplying by Q we have /(2;)(5(2;) = P{z). Now equating the coefficients on both sides 
we get the following system of linear equations: 

mo^o = Aq 

mo-Bi + miSo = Ai 



uioBn-i + miBN-2 H h tun-iBq = A^^i 

mo-Biv + t^iBn-i H h niN-iBi + mivi?o = 

miSiv + rn2BM-i H h ttinBi + tiin+iBq = 



The rest of the equations in this system are obtained by further shifts of the indices of the 
moments, and so they form a Hankel-type matrix. 

Now, being a rational function of degree A^, I{z) is uniquely defined by its first 2N Taylor 
coefficients (the difference of two such functions cannot vanish at zero with the order higher 
than 2N — 1). We conclude that the linear system consisting of the first 2N equations as 
above is uniquely solvable up to a common scaling of P and Q (of course, this fact follows 
also form a general Pade approximation theory - see |22j). 

Now a solution procedure for the Prony system can be described as follows: 

1. Solve a linear system of the first 2N equations as above (with the coefficients - the 
known moments m^) to find the moments generating function I{t) in the form I{z) = ^||y. 

2. Represent I{z) in a standard way as the sum of elementary fractions I{z) = XljLi i^fc- 
(Equivalently, find poles and residues of I{z)). Besides algebraic operations, this requires 
just finding the roots of the polynomial Q{z). Then (a^, x^), j = 1, . . . , N form the unique 
solution of the Prony system (fT2!) . 

3.2. Multi-dimensional case. Let us recall our multi-dimensional notations: x = (xi, . . . , Xd) G 
M'^, j is a scalar index, while k = {ki, . . . ,kd), i = {ii, . . . ,id) and n = (ni, . . . , Ud) are 
multi-indices. Partial ordering of multi- indices is given hj k < k' kp < k'^, p = 1, . . . ,d. 

Let x^ = (x{, . . . , x^l|), j = 1, . . . , N. With the above notation, the multidimensional Prony 
system has exactly the same form as in the one-dimensional case: 

N 

(14) J2 %(^^)^ = ^k, e A; > (0, . . . , 0). 

i=i 

Exactly as above we get that for z = {zi, . . . , Zd) E the moments generating function 
H^) = XlfceNd 'iT^kZ^ is a rational function of degree Nd of the form 
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N d 



1 



(15) 



j=l 1=1 



Representing I{z) as we get exactly in the same way as above an infinite system of 
linear equations for the coefficients of P and Q, with a Hankel-type matrix formed by the 
moments m^. By the same consideration as above, after we take enough equations in this 
system the solution is unique up to a rescaling (see [H |22] and references therein). 

However, from this point the multi-dimensional situation becomes essentially more com- 
plicated. While in dimension one I{z) can be, essentially, any rational function of degree 
(naturally represented as the sum of elementary fractions), in several variables I{z) has 
a very special form given by f|T5|) . This fact can be easily understood via counting the de- 
grees of freedom: the number of unknowns in f|T^ is N{d + 1) while a rational function 
of degree Nd in d variables has much more coefficients than that, for d > 1. It would be 
desirable to use as many equations from f lT^ as the number of unknowns, but the method 
outlined above ignores a specific structure of I{z) and requires as many equations as in a 
Fade reconstruction for a general rational function of degree Nd. 

We treat this problem in [28] analyzing the structure of singularities of I{z) and on this 
base proposing a robust reconstruction algorithm. Let us give here one simple special case 
of this algorithm, which separates the variables in the problem. 

3.2.1. Separation of variables in the multi- dimensional Prony system. Let us assume that 
we know a priori that the solution {aj,x^), j = 1, . . . , N of multi-dimensional Prony system 
(fT4|) is such that all the coordinates ocj, I = l,...,d of the points x\ j = 1,...N are 
pairwise distinct. Moreover, we assume that 7^ aj^ for ji 7^ j2- Under these assumptions 
we proceed as follows: 

Consider "partial moment generating functions" Im(t), t G C, m = 1, . . . , d, defined by 



where Cm is a multi-index with {em)j = for m ^ j and 1 otherwise. We have the following 
simple fact: 

Proposition 2. Im{t) is a one- dimensional moments generating function of the form 



00 



(16) 




r=l 



(17) 




It coincides with the restriction of I{z) to the m-th coordinate axis in C^. 
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Proof. Let us evaluate I{z) along the m-th coordinate axis, that is on the line z = tcm with 
Cm as above and t G C. We get 



N d ^ ^1 



j = l 1 = 1 ^ -^l ^\^m)l 

which is the moments generating function ( ITTl) . Now, to express Imit) through the multi- 
dimensional moments ruk we notice that any monomial x'^, with k ^ pcm vanishes identically 
on the m-th coordinate axis. Hence 



r=0 



This shows that Imit) = I{tem) and completes the proof of the proposition. □ 

Now applying the method described in Section fl3.ip above we find for each m = 1, . . . , d 
the coordinates x^, . . . , x^ and (repeatedly) the coefficients ai, . . . , oat. It remains to arrange 
these coordinates into the points ^). This presents a certain combinatorial 

problem, since Prony system f|T^ is invariant under permutations of the index j. Under 
the assumptions above we proceed as follows: for each m = 1, . . . , d we have obtained the 
(unordered) collection of the pairs (a^, x^), j = ^, ■ ■ ■ , N. By assumptions a^^ ^ aj^ for ji ^ 
j2- Hence we can arrange in a unique way all the pairs (a^, x^), i = 1, ■ • • , m = 1, . . . , 
into the string [(ai,xj), . . . , (ai,x^)], . . . , [(aAr,xf^), . . . , (aAr,x^)]. This gives us the desired 
solution of multi-dimensional Prony system ([1] 



Notice that the assumption a^^ ^ aj^ for ji ^ j2 is essential here. Indeed, for xi 7^ X2 
and x^ = (xi,X2), x^ = (x2,xi), = (xi,xi), x^ = (x2,X2) we have = (x"^)'^ -f (x^)*^ = 
^^i-jfc _l_ (■^2^fc ^ each of the coordinate axes. So the (unique up to permutations of the 
index j) solution of the Prony system cannot be reconstructed from these moments only. 

Another remark is that the separation of variables as described above requires knowledge 
of 2dN moments rrik {2N on each of the coordinate axes). This is almost twice more than 
N{d + 1) unknowns. This number can be significantly reduced in some cases. See |28] for 
further investigation in both of these directions. 

Stability estimates for the solution of one-dimensional Prony system (Section Ej) can be 
easily extended to the case considered in the present subsection. We do not give here explicit 
statement of this result. 



4. Reconstruction of piecewise D-finite functions from moments 

In this section we present an overview of our previous findings on algebraic reconstruction 
of a certain general class of signals. See [3] for the complete details. 

Let g : [a,b] ^ M. consist of -f- 1 "pieces" go, ■ ■ ■ gN with A > jump points 



a = Xo < Xi . . . < Xat < X^r+i = b 
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Furthermore, let g satisfy on each continuity interval some linear homogeneous differential 
equation with polynomial coefficients: D gn = 0, n = 0, . . . , N where 

(18) ® = E(E^m^O^ (^^.^^) 

Each gn may be therefore written as a linear combination of functions which are a 

basis for the space = {/ : D / = 0}: 

r 

(19) gn{x) = '^ai^nUi{x), n = 0,l,...,N 

1=1 

We term such functions g "piecewise D-finite" . Many "simple" functions may be represented 
in this framework, such as polynomials, trigonometric polynomials and algebraic functions. 
The sequence {m^ = mi^{g)} is given by the usual moments 

mk{g) = / x''g{x)dx 

J a 

Piecewise D-finite Reconstruction Problem. Given r, {fcj}, A^, a, b and the moment se- 
quence {rrik} of a piecewise D -finite function g, reconstruct all the parameters {bij}, {xi}, {ai,n}- 

The above problem can be solved as follows. 

(1) It turns out that the moment sequence of every piecewise D- finite function g satisfies 
a linear recurrence relation, such that the coefficients of D annihilating every piece of 
g and the discontinuity locations {xj} can be recovered from the moments by solving 
a linear systems of equations plus a nonlinear step of polynomial roots finding. As 
we shall explain below, in many cases this step is equivalent to solving a certain 
generalized form of the previously mentioned Prony system f|T2|) . 

(2) The function g itself can be subsequently reconstructed by numerically calculating 
the basis for the space A/j) and solving an additional linear system of equations to 
recover the coefficients 

The above algorithm has been tested on reconstruction of piecewise polynomials, piecewise 
sinusoids and rational functions - see [3] for details. The results reported there were relatively 
accurate for low noise levels. In this paper we continue to explore the numerical stability of 
the method - see Sections |5] and [61 

Now let us show how the Prony system arises in the piecewise D-finite reconstruction 
method. Consider the distribution 2) g. It is of the form (see [Sj Theorem 2.12]) 



N 

(20) S)^7 = 5^5^a,„,5« 

j=l i=0 

where Xj G M are the discontinuity locations, G M are the associated "jump magnitudes" 
which depend on the values {g''^*^(x^)}, and 5{x) is the Dirac 5-function. In particular, when 
g is just a piecewise constant function, we have 2) = ^ and so 

N 

(21) D g = '^aj5{x - Xj) 
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where aj = g{x'j') — g{xj ). Let us now apply the moment transform to the equations (l2Ti) 
and ( 12U]) . We get, correspondingly, 

N 

(22) mfc(2)^) = J^a^x^^; 

(23) mfc(D^) = ^ ^ a,jA;(A; - 1) X ■ ■ ■ X (A; - i + l)x^^"\ 

j=l i=Q 

The system fl22p is of course identical to the previously considered ( I12p . However, the 
following question arises: how are the numbers in the left-hand side of fl22p and f l23p related 
to the known quantities nik ((/)? It turns out that the numbers ruk (2)^?) are certain linear 
combinations of these known moments, with coefficients which are determined by 2) in a well- 
defined way. The conclusion is as follows: if the operator T) which annihilates every piece 
of g is known (but the other parameters are no^), then one can recover the discontinuity 
locations of g by solving the Prony-like system fl23l) . We term the system fl23l) "confiuent 
Prony system". It can be solved in a similar manner to the standard Prony system. A 
unique solution exists whenever all the Xj S QjTQ distinct and the highest coefficients a; 
are nonzero. We provide the details in [2]. 

5. Accuracy of solving Prony-like systems 

In the preceding sections we have shown that several algebraic methods for nonlinear 
reconstruction can be reduced to solving certain types of systems of equations, the most basic 
one of which is the Prony system (fT2|) . (122|) . A crucial factor for the eventual applicability of 
the reconstruction methods is the sensitivity of solving these systems to measurement errors. 

In this section we consider two such systems - (122|) and (!23|) and provide some theoretical 
results regarding the local sensitivity of their solution to noise. These results will hopefully 
enable further "global" analysis. 

We consider the following question: if the measurements in the left-hand side of fl22p are 
known with error at most e, how well can one recover the unknowns {aj,Xj}7 We regard this 
stabihty problem to be absolutely central in Algebraic Sampling. To our best knowledge, 
no general treatment of this problem exists, therefore we consider our results below to be a 
step in this direction. 

For simplicity, let us assume that the number of equations in fl22|) equals the number of 
unknowns, which in this case is P = 2A^. Let us further consider the "measurement map" 
V -.R^ given by ([22]) (we call it the "Prony map"): 

V{{aj},{xj}) = {nik}^-^ 

Then, one possible answer (by no means a complete one) to the question posed above can 
be given in terms of the local Lipschitz constant of the "solution map" V~^, whenever this 
inverse is defined. We then obtain the following result. 



This assumption is realistic, for instance when reconstructing piecewise-polynomials or piecewise- 
sinusoids. 
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Theorem 3. Let {mk}k=o,...,p-i be the exact unperturbed moments of the model (12T]) . 74s- 
sume that all the xj 's are distinct and also aj ^ for j = 1, . . . , N. Now let irfk be pertur- 
bations of the above moments such that maxfe \mk — mk\ < e. Then, for sufficiently small e, 
the perturbed Prony system has a unique solution which satisfies: 

\Xj "^jl — ^1 ^ I ^ J I 
[o-j — fljl < CiE 

where Ci is an explicit constant depending only on the geometry of xi, . . . ,xn. 

Proof. Consider the Jacobian determinant of the map V. It is easily factorized as foUows: 





■ 1 





1 


" 






dV = 


Xi 


1 


Xn 


1 


X diag{Di, . 






p 

_x 


Fxf-i ... 


p 




✓ 








dcf, 

= y(xi,. 


.,xjv) 









where is a 2 x 2 block 



def 



1 

tti 



The matrix = y (xi, . . . , xat) is a special case of the so-called confluent Vandermonde 
matrix, well-known in numerical analysis ([HI [111 [131 [E])- In particular, the paper 
Theorem 3] contains the following estimate for the norm of V^^: 



\V-^\\oo < max h, W 

l<i<N J-J- \\Xi — 



n 




1 -^i 





where 



, def 

0,- = max 



(■ 



+ IxJ, 1 + 2(1 



■i)E 



Xi 



Now since the x^-'s are distinct and aj ^ 0, the Jacobian is non-singular and so in a sufficiently 
small neighbourhood of the exact solution, the map V is approximately linear. By the inverse 
function theorem 

dp-^ 



and so taking C\ =^ ||y ^ 



diag{Dr\...,D^^}xF-i 
completes the proof. 



□ 



By a similar technique with slightly more involved calculations we obtain the following 
result for the system ( 123|) . 

Theorem 4. Let {mfc}fc=o,...,p-i be the exact unperturbed moments of the model fl20|) . where 
P = XljLi + ^- Assume that all the Xj 's are distinct and also ai.-ij ^ for j = 1, . . . , N . 
Now let ink be perturbations of the above moments such that max^ \mk — rrTkl < Then, 
for sufficiently small e, the perturbed confluent Prony system has a unique solution which 
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satisfies: 



I ^« ,i 




l<i<lj-l 



where C2 is a constant depending only on the nodes xi, . . . ,xn and the multiplicities li, . . . ,1^- 

Proof outline. As before, we obtain a factorization of the Jacobian determinant dP as a 
product of a confluent Vandermonde matrix V{xi, li, . . . ,xn,In) (defined in a similar manner 
but with each column having Ij — 2 "confluencies" ) and the block diagonal matrix D = 
diagl^i, . . . , -Dat} where 

"1 ••• 

1 ■■■ ao,: 

^3 

_0 ■ ■ ■ cti^-ij. 

Inverting this expression and taking C2 = ||^~^||oo completes the proof. □ 
An important generalization would be to consider the mappings 

({a,,},{a;,})"=^{m,}r=r' 

and investigate the reconstruction error as s — i- 00. Such a formulation makes sense in the 
Fourier reconstruction problem ( |H [9l [H] ) . This is a work in progress and we plan to present 
the results separately ([2])- 

6. Numerical experiments 

We have tested the piecewise D-finite reconstruction method on a simple case of a piecewise- 
constant signal. The implementation details are identical to those used in jSj Appendix] . As 
can be seen from Figured! the reconstruction is accurate even in the presence of medium- level 
noise. 

We have already mentioned that for a piecewise constant function /, the distribution /' is 
of the form fl2T]) . Therefore, the piecewise D-finite reconstruction problem for / essentially 
reduces to solving the Prony system ( l22l) . and so the local estimates of Theorem [3] should 
apply in the case of a very small noise level. This prediction is partially confirmed by the 
results of our second experiment, presented in Figure O In particular, it can be seen that 
indeed |Axj| ~ while the accuracy of other jump points \^Xi\ does not depend on 

\aj\ for j 7^ i. While this is certainly an encouraging result, more investigation is clearly 
needed in order to fully understand the dependence of the error on all the parameters of the 
problem. Such an investigation should, in our opinion, concentrate on the global structure 
of the mapping V. 
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